R Notebook

Dependencies

# data wrangling
library(tidyverse)

# spatial data wrangling
library(sf)

# data visualisation
library(viridis) 

# format data visualisations
library(ggthemes)
library(patchwork)
library(showtext)
library(scales)
library(classInt)
library(ggtext)

# create maps
library(leaflet)
library(tmap)
library(mapdeck)
library(patchwork)
library(cowplot)

Setting theme

Set font style

# clean workspace
rm(list=ls())
# load font
font_add_google("Roboto Condensed", "robotocondensed")
# automatically use showtext to render text
showtext_auto()

Theme for maps

theme_map <- function(...) {
  theme_tufte() +
  theme(
    text = element_text(family = "robotocondensed", size = 20),
    # remove all axes
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank()
    )
}

Theme for plots

theme_tufte2 <- function(...) {
  theme_tufte() +
  theme(
    text = element_text(family = "robotocondensed", size = 20),
    )
}

Argentina

Data

mobility data

## need to read this in 4 separate times - 202

# 2020 out
df20_b <- readRDS("/Volumes/RECAST/data/outputs/argentina/movements/2020_04_mov.rds") %>% 
  mutate(GEOMETRY = NULL) %>% 
  dplyr::filter(country == "AR") %>% 
  st_as_sf(coords = c("end_lon", "end_lat"), 
                                      crs = 'EPSG:4326')

# 2022 out
df22_b <- readRDS("/Volumes/RECAST/data/outputs/argentina/movements/2022_03_mov.rds") %>% 
  mutate(GEOMETRY = NULL) %>% 
  dplyr::filter(country == "AR") %>% 
  st_as_sf(coords = c("end_lon", "end_lat"), 
                                      crs = 'EPSG:4326')

boundary data

adm_shp <- st_read("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_ARG_shp/gadm41_ARG_2.shp") %>% 
  st_simplify(preserveTopology = T,
              dTolerance = 1000) %>% 
  st_make_valid() %>% 
  st_transform(crs = 'EPSG:4326')

Reading layer `gadm41_ARG_2' from data source 
  `/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_ARG_shp/gadm41_ARG_2.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 502 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -73.56056 ymin: -55.06153 xmax: -53.59184 ymax: -21.78137
Geodetic CRS:  WGS 84

adm1_shp <- st_read("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_ARG_shp/gadm41_ARG_1.shp") %>% 
  st_simplify(preserveTopology = T,
              dTolerance = 1000) %>% 
  st_make_valid() %>% 
  st_transform(crs = 'EPSG:4326')

Reading layer `gadm41_ARG_1' from data source 
  `/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_ARG_shp/gadm41_ARG_1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 24 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -73.56056 ymin: -55.06153 xmax: -53.59184 ymax: -21.78137
Geodetic CRS:  WGS 84

region_shp <- st_read("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_ARG_shp/gadm41_ARG_0.shp")%>% 
  st_simplify(preserveTopology = T,
              dTolerance = 1000) %>% 
  st_make_valid() %>% 
  st_transform(crs = 'EPSG:4326')

Reading layer `gadm41_ARG_0' from data source 
  `/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_ARG_shp/gadm41_ARG_0.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -73.56056 ymin: -55.06153 xmax: -53.59184 ymax: -21.78137
Geodetic CRS:  WGS 84

Data wrangling

area_names <- c("José C. Paz", "San Miguel", "Morón", "Capital Federal", "Esteban Echeverría", "Florencio Varela")

# filter for the area names which are within Buenos Aires
df20_out <- df20_b %>% dplyr::filter(start_polygon_name %in% area_names)
df22_out <- df22_b %>% dplyr::filter(start_polygon_name %in% area_names)

df20_out$start_polygon_name <- 'Buenos Aires'
df22_out$start_polygon_name <- 'Buenos Aires'

df20_out <- df20_out %>% filter(!end_polygon_name %in% area_names)
df22_out <- df22_out %>% filter(!end_polygon_name %in% area_names)

Define distance

df20_out <- df20_out %>% mutate(
  distance_class = case_when(length_km < 100 ~ "<100",
                             length_km >= 100 ~ ">100"))

df22_out <- df22_out %>% mutate(
  distance_class = case_when(length_km < 100 ~ "<100",
                             length_km >= 100 ~ ">100"))

preserve a unique geometry

df20_out <- df20_out %>%
  group_by(end_polygon_name) %>%
  mutate(geometry = st_union(geometry)) %>%
  ungroup()

df22_out <- df22_out %>%
  group_by(end_polygon_name) %>%
  mutate(geometry = st_union(geometry)) %>%
  ungroup()

Sum of flows

outflows_df_20 <- df20_out %>% 
  filter(start_polygon_name != end_polygon_name) %>% 
  group_by(start_polygon_name, end_polygon_name, geometry, distance_class) %>% 
  dplyr::summarise(
    sum_outflow = sum(n_crisis, na.rm = T)) %>% 
  ungroup() 

`summarise()` has grouped output by 'start_polygon_name', 'end_polygon_name', 'geometry'. You can override using the `.groups` argument.

outflows_df_22 <- df22_out %>% 
  filter(start_polygon_name != end_polygon_name) %>% 
  group_by(start_polygon_name, end_polygon_name, geometry, distance_class) %>% 
  dplyr::summarise(
    sum_outflow = sum(n_crisis, na.rm = T)) %>% 
  ungroup() 

`summarise()` has grouped output by 'start_polygon_name', 'end_polygon_name', 'geometry'. You can override using the `.groups` argument.

Check alignment between data frames

p <- ggplot() + 
  geom_sf(data = adm_shp,
          color = "gray60", 
          size = 0.1) 

last_plot()

p <- p +
  geom_point(data = outflows_df_20,
    aes(geometry = geometry),
    stat = "sf_coordinates"
  ) 
last_plot()

Add year

outflows_df_20$year <- '2020'
outflows_df_22$year <- '2022'

Bind outflows

outflows_df <- rbind(outflows_df_20, outflows_df_22)

Join flows to polygons

mob_indicators_1 <- st_join(adm_shp, outflows_df)

classify data into quantiles

mob_indicators_1$new_jenk_class <- classify_intervals(mob_indicators_1$sum_outflow, n = 5, style = "quantile", factor = TRUE)

Warning: var has missing values, omitted in finding classes

outflow_labels_1 <- levels(mob_indicators_1$new_jenk_class)
outflow_labels_1 <- gsub("^.|.$", "",  outflow_labels_1)
outflow_labels_1 <- gsub("\\.[0-9]+", "",  outflow_labels_1)
outflow_labels_1 <- gsub("\\,", "-",  outflow_labels_1)
levels(mob_indicators_1$new_jenk_class) <- outflow_labels_1

# mob_indicators_1 <- mob_indicators_1 %>% 
  # mutate(new_jenk_class = str_replace(new_jenk_class, ",", "-"))

# change geometry
shp_reg <- region_shp %>% st_transform(crs = 'EPSG:4326')

mob_indicators_1$new_jenk_class <- mob_indicators_1$new_jenk_class %>% replace_na("10-32")

mob_indicators_1 <- na.omit(mob_indicators_1)

Plotting

Creating individual country map

shp_reg$centroid <- shp_reg %>% 
  st_centroid() %>% 
  st_geometry()

Warning: st_centroid assumes attributes are constant over geometries

padding_width <- 17 
padding_height <- 22

ggplot(data = mob_indicators_1, 
                           aes(fill = new_jenk_class )) +
  geom_sf(col = "white", size = .2) + 
  coord_sf() +
  scale_fill_brewer(palette = "OrRd", 
                    direction = 1, 
                    labels = c("10-32", "32-242", "243-1,640", "1641-30,291", "> 30,292"), 
                    na.value="black") +
  facet_grid(distance_class ~ year) +
  labs(title = "A. Argentina",
       fill = "Number of out-moves") +
  theme_map() +
  guides(
    color='none',
    fill = guide_legend(
      keywidth = 4, 
      keyheight = 1,
      nrow = 1,
      title.position="top",
      label.position="bottom"
    )
  ) +
  theme(plot.title = element_text(size = 16, face = "bold"),
        plot.margin=margin(1,0,1,0,"cm"),
        legend.title = element_markdown(
          size=10, face = "bold", hjust=0.5, lineheight=0.45,
          color="black",
          margin=margin(0,0,-0.2,0,"cm")
          ),
        legend.text = element_text(size = 8),
        legend.position = "bottom",
        legend.spacing.x = unit(0, 'cm'),
        panel.background = element_rect(fill = "gray98", colour = "gray98")
        ) +
    geom_sf(data = shp_reg,
          col = "black", 
          size = 1,
          fill = "transparent") +
            coord_sf(xlim = c(shp_reg$centroid[[1]][1] - padding_width, 
                              shp_reg$centroid[[1]][1] + padding_width), 
           #          ylim = c(shp_reg$centroid[[1]][2] - padding_height, 
           #                   shp_reg$centroid[[1]][2] + padding_height), 
                    expand = FALSE)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

outflow_plot_arg <- last_plot()

Saving map

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/argentina/cloropleth-arg.png", units="in", width=8, height=10, res=300)
  outflow_plot_arg
dev.off()

null device 
          1 

Create map for Buenos Aires

Short-distance

ba_deparment_shp <- adm1_shp %>% filter(NAME_1 == "Buenos Aires" |
                                          NAME_1 == "Ciudad de Buenos Aires" |
                                          NAME_1 == "Entre Ríos")

ba_deparment_shp$centroid <- ba_deparment_shp %>% 
  st_centroid() %>% 
  st_geometry()

Warning: st_centroid assumes attributes are constant over geometries

padding_width <- 3.5 
padding_height <- 3.2

# 2020
mob_indicators_1 %>% filter(NAME_1 == "Buenos Aires"| 
                              NAME_1 == "Ciudad de Buenos Aires" | 
                              NAME_1 == "Entre Ríos" |
                              NAME_1 == "Córdoba" |
                              NAME_1 == "Santa Fe") %>% 
  filter(year == "2020") %>% 
  filter(distance_class == "<100") %>% 
   ggplot(aes(fill = new_jenk_class)) +
     geom_sf(col = "white", size = .1) +
   coord_sf() +
  scale_fill_brewer(palette = "OrRd", 
                    direction = 1, 
                    labels = c("10-32", "32-242", "243-1,640", "1641-30,291", "> 30,292")) +
   theme_map() +
   theme(plot.title = element_text(size = 60),
         legend.position = "none",
         panel.background=element_rect(colour="black")
         ) +
   labs(title = "Buenos Aires") +
    geom_sf(data = ba_deparment_shp,
          col = "black", 
          size = 1,
          fill = "transparent") +
            coord_sf(xlim = c(ba_deparment_shp$centroid[[1]][1] - -.5, 
                              ba_deparment_shp$centroid[[1]][1] + padding_width), 
                     ylim = c(ba_deparment_shp$centroid[[1]][2] - -.5, 
                              ba_deparment_shp$centroid[[1]][2] + padding_height), 
                    expand = FALSE)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

ba_map2020 <- last_plot()

# 2022
mob_indicators_1 %>% filter(NAME_1 == "Buenos Aires"| 
                              NAME_1 == "Ciudad de Buenos Aires" | 
                              NAME_1 == "Entre Ríos" |
                              NAME_1 == "Córdoba" |
                              NAME_1 == "Santa Fe") %>% 
  filter(year == "2022") %>% 
  filter(distance_class == "<100") %>% 
   ggplot(aes(fill = new_jenk_class)) +
     geom_sf(col = "white", size = .1) +
   coord_sf() +
  scale_fill_brewer(palette = "OrRd", 
                    direction = 1, 
                    labels = c("10-32", "32-242", "243-1,640", "1641-30,291", "> 30,292")) +
   theme_map() +
   theme(plot.title = element_text(size = 60),
         legend.position = "none",
         panel.background=element_rect(colour="black")
         ) +
   labs(title = "Buenos Aires") +
    geom_sf(data = ba_deparment_shp,
          col = "black", 
          size = 1,
          fill = "transparent") +
            coord_sf(xlim = c(ba_deparment_shp$centroid[[1]][1] - -.5, 
                              ba_deparment_shp$centroid[[1]][1] + padding_width), 
                     ylim = c(ba_deparment_shp$centroid[[1]][2] - -.5, 
                              ba_deparment_shp$centroid[[1]][2] + padding_height), 
                    expand = FALSE)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

ba_map2022 <- last_plot()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/argentina/cloropleth-ba2020.png", units="in", width=8, height=10, res=300)
  ba_map2020
dev.off()

null device 
          1 

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/argentina/cloropleth-ba2022.png", units="in", width=8, height=10, res=300)

  ba_map2022
dev.off()

null device 
          1 

Create bar plot

Short-distance

col_pal <- c(#"#fef0d9",
             #"#fdcc8a",
             #"#fc8d59",
             "#e34a33",
             "#b30000")

# 2020
mob_indicators_1 %>% filter(year == "2020") %>%
  filter(distance_class == "<100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
    distinct(end_polygon_name, .keep_all=TRUE) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(end_polygon_name, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
       # axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_short2020 <- last_plot()

# 2022
mob_indicators_1 %>% filter(year == "2022") %>%
  filter(distance_class == "<100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  distinct(end_polygon_name, .keep_all=TRUE) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(end_polygon_name, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
        #axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_short2022 <- last_plot()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/argentina/barp_short2020.png", units="in", width=8, height=10, res=300)
  barp_short2020
dev.off()

null device 
          1 

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/argentina/barp_short2022.png", units="in", width=8, height=10, res=300)

  barp_short2022
dev.off()

null device 
          1 

Long-distance

2020

col_pal <- c("#fef0d9",
             "#fdcc8a"
             #"#fc8d59"
             #"#e34a33",
             #"#b30000"
  )

mob_indicators_1 %>% filter(year == "2020") %>%
  filter(distance_class == ">100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(NAME_2, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
       # axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_long2020 <- last_plot()

2022

col_pal <- c(#"#fef0d9",
             #"#fdcc8a"
             "#fc8d59",
             "#e34a33"
             #"#b30000"
  )

mob_indicators_1 %>% filter(year == "2022") %>%
  filter(distance_class == ">100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(NAME_2, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
       # axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_long2022 <- last_plot()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/argentina/barp_long2020.png", units="in", width=8, height=10, res=300)
  barp_long2020
dev.off()

null device 
          1 

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/argentina/barp_long2022.png", units="in", width=8, height=10, res=300)

  barp_long2022
dev.off()

null device 
          1 

Chile

Data

mobility data

## need to read this in 4 separate times - 202

# 2020 out
df20_b <- readRDS("/Volumes/RECAST/data/outputs/chile/movements/2020_04_mov.rds") %>% 
  mutate(GEOMETRY = NULL) %>% 
  dplyr::filter(country == "CL") %>% 
  st_as_sf(coords = c("end_lon", "end_lat"), 
                                      crs = 'EPSG:4326')

# 2022 out
df22_b <- readRDS("/Volumes/RECAST/data/outputs/chile/movements/2022_03_mov.rds") %>% 
  mutate(GEOMETRY = NULL) %>% 
  dplyr::filter(country == "CL") %>% 
  st_as_sf(coords = c("end_lon", "end_lat"), 
                                      crs = 'EPSG:4326')

boundary data

adm_shp <- st_read("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/Chile_shp/adm/province/PROVINCIAS_2020.shp")%>% 
  st_simplify(preserveTopology = T,
              dTolerance = 1000) %>% 
  st_make_valid() %>% 
  st_transform(crs = 'EPSG:4326')

Reading layer `PROVINCIAS_2020' from data source 
  `/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/Chile_shp/adm/province/PROVINCIAS_2020.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 56 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -109.4488 ymin: -56.51273 xmax: -66.42812 ymax: -17.4984
Geodetic CRS:  SIRGAS-Chile 2002

adm1_shp <- st_read("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/Chile_shp/adm/region/REGIONES_2020.shp")%>% 
  st_simplify(preserveTopology = T,
              dTolerance = 1000) %>% 
  st_make_valid() %>% 
  st_transform(crs = 'EPSG:4326')

Reading layer `REGIONES_2020' from data source 
  `/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/Chile_shp/adm/region/REGIONES_2020.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 16 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -109.4488 ymin: -56.51273 xmax: -66.42812 ymax: -17.4984
Geodetic CRS:  SIRGAS-Chile 2002

region_shp <- st_read("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/Chile_shp/adm/country/gadm41_CHL_0.shp")%>% 
  st_simplify(preserveTopology = T,
              dTolerance = 1000) %>% 
  st_make_valid() %>% 
  st_transform(crs = 'EPSG:4326')

Reading layer `gadm41_CHL_0' from data source 
  `/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/Chile_shp/adm/country/gadm41_CHL_0.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -109.4549 ymin: -55.98 xmax: -66.41821 ymax: -17.49859
Geodetic CRS:  WGS 84

Data wrangling

df20_out <- df20_b %>% 
  filter(start_polygon_name == 'Santiago')
  
df22_out <- df22_b %>% 
  filter(start_polygon_name == 'Santiago')

Define distance

df20_out <- df20_out %>% mutate(
  distance_class = case_when(length_km < 100 ~ "<100",
                             length_km >= 100 ~ ">100"))

df22_out <- df22_out %>% mutate(
  distance_class = case_when(length_km < 100 ~ "<100",
                             length_km >= 100 ~ ">100"))

preserve a unique geometry

df20_out <- df20_out %>%
  group_by(end_polygon_name) %>%
  mutate(geometry = st_union(geometry)) %>%
  ungroup()

df22_out <- df22_out %>%
  group_by(end_polygon_name) %>%
  mutate(geometry = st_union(geometry)) %>%
  ungroup()

Sum of flows

outflows_df_20 <- df20_out %>% 
  filter(start_polygon_name != end_polygon_name) %>% 
  group_by(start_polygon_name, end_polygon_name, geometry, distance_class) %>% 
  dplyr::summarise(
    sum_outflow = sum(n_crisis, na.rm = T)) %>% 
  ungroup() 

`summarise()` has grouped output by 'start_polygon_name', 'end_polygon_name', 'geometry'. You can override using the `.groups` argument.

outflows_df_22 <- df22_out %>% 
  filter(start_polygon_name != end_polygon_name) %>% 
  group_by(start_polygon_name, end_polygon_name, geometry, distance_class) %>% 
  dplyr::summarise(
    sum_outflow = sum(n_crisis, na.rm = T)) %>% 
  ungroup() 

`summarise()` has grouped output by 'start_polygon_name', 'end_polygon_name', 'geometry'. You can override using the `.groups` argument.

Check alignment between data frames

p <- ggplot() + 
  geom_sf(data = adm_shp,
          color = "gray60", 
          size = 0.1) 

last_plot()

p <- p +
  geom_point(data = outflows_df_20,
    aes(geometry = geometry),
    stat = "sf_coordinates"
  ) 
last_plot()

Add year

outflows_df_20$year <- '2020'
outflows_df_22$year <- '2022'

Bind outflows

outflows_df <- rbind(outflows_df_20, outflows_df_22)

Join flows to polygons

mob_indicators_1 <- st_join(adm_shp, outflows_df)

classify data into quantiles

mob_indicators_1$new_jenk_class <- classify_intervals(mob_indicators_1$sum_outflow, n = 5, style = "quantile", factor = TRUE)

Warning: var has missing values, omitted in finding classes

outflow_labels_1 <- levels(mob_indicators_1$new_jenk_class)
outflow_labels_1 <- gsub("^.|.$", "",  outflow_labels_1)
outflow_labels_1 <- gsub("\\.[0-9]+", "",  outflow_labels_1)
outflow_labels_1 <- gsub("\\,", "-",  outflow_labels_1)
levels(mob_indicators_1$new_jenk_class) <- outflow_labels_1

# mob_indicators_1 <- mob_indicators_1 %>% 
  # mutate(new_jenk_class = str_replace(new_jenk_class, ",", "-"))

# change geometry
shp_reg <- region_shp %>% st_transform(crs = 'EPSG:4326')

mob_indicators_1$new_jenk_class <- mob_indicators_1$new_jenk_class %>% replace_na("10-420")

mob_indicators_1 <- na.omit(mob_indicators_1)

Plotting

Crop bounding box

bbox_new <- st_bbox(adm_shp) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

bbox_new[1] <- bbox_new[1] + (0.6 * xrange) # xmin - left

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

ggplot() + 
  geom_sf(data = adm_shp,
          color = "gray60", 
          size = 0.1) +
  geom_point(data = outflows_df_20,
    aes(geometry = geometry),
    stat = "sf_coordinates",
    size = .1
  ) +
  coord_sf(xlim = st_coordinates(bbox_new)[c(1,2),1], # min & max of x values
           ylim = st_coordinates(bbox_new)[c(2,3),2]) + # min & max of y values
  theme_void()

Creating individual country map

shp_reg$centroid <- shp_reg %>% 
  st_centroid() %>% 
  st_geometry()

Warning: st_centroid assumes attributes are constant over geometries

padding_width <- 17
padding_height <- 22

ggplot(data = mob_indicators_1, 
                           aes(fill = new_jenk_class )) +
  geom_sf(col = "transparent", size = .2) + 
  coord_sf() +
  scale_fill_brewer(palette = "OrRd", 
                    direction = 1, 
                    labels = c("10-421", "422-1231", "1,232-3,370", "3,371-22,761", "> 22,761"),
                    na.value="black") +
  facet_grid(distance_class ~ year) +
  labs(title = "B. Chile",
       fill = "Number of out-moves") +
  theme_map() +
  guides(
    color='none',
    fill=guide_legend(
      keywidth = 4,
      keyheight = 1,
      nrow = 1,
      title.position="top",
      label.position="bottom"
    )
  ) +
  theme(plot.title = element_text(size = 16, face = "bold"),
        plot.margin=margin(1,0,1,0,"cm"),
        legend.title = element_markdown(
          size=10, face = "bold", hjust=0.5, lineheight=0.45,
          color="black",
          margin=margin(0,0,-0.2,0,"cm")
          ),
        legend.text = element_text(size = 9),
        legend.position = "bottom",
        legend.spacing.x = unit(0, 'cm'),
        panel.background = element_rect(fill = "gray99", colour = "gray99")
        ) +
    geom_sf(data = region_shp,
          col = "black", 
          size = 1,
          fill = "transparent") +
            coord_sf(xlim = c(shp_reg$centroid[[1]][1] - padding_width, 
                              shp_reg$centroid[[1]][1] + padding_width), 
                     ylim = c(shp_reg$centroid[[1]][2] - padding_height, 
                              shp_reg$centroid[[1]][2] + padding_height), 
                    expand = FALSE)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

outflow_plot_chi <- last_plot()

Saving country map

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/chile/cloropleth-chi.png", units="in", width=8, height=10, res=300)
  last_plot()
dev.off()

null device 
          1 

Create map for Santiago

Short-distance

metro_region_shp <- adm1_shp %>% filter(REGION == "Metropolitana de Santiago" |
                                         REGION == "Valparaíso" | 
                                         REGION == "Libertador General Bernardo O'Higgins")

metro_region_shp$centroid <- metro_region_shp  %>% 
  st_centroid() %>% 
  st_geometry()

Warning: st_centroid assumes attributes are constant over geometries

padding_width <- 2
padding_height <- 1.8

col_pal <- c(#"#fef0d9",
             #"#fdcc8a",
             "#fc8d59",
             "#e34a33",
             "#b30000")

# 2020
mob_indicators_1 %>% filter(REGION == "Metropolitana de Santiago" |
                                         REGION == "Valparaíso" | 
                                         REGION == "Libertador General Bernardo O'Higgins") %>% 
  filter(year == "2020") %>% 
  filter(distance_class == "<100") %>% 
   ggplot(aes(fill = new_jenk_class)) +
     geom_sf(col = "white", size = .1) +
   coord_sf() +
   scale_fill_manual(values= col_pal) +
  # scale_fill_brewer(palette = "OrRd", 
  #                   direction = 1, 
  #                   labels = c("10-421", "422-1231", "1,232-3,370", "3,371-22,761", "> 22,761")) +
   theme_map() +
   theme(plot.title = element_text(size = 60),
         legend.position = "none",
         panel.background=element_rect(colour="black")
         ) +
   labs(title = "Santiago") +
    geom_sf(data = metro_region_shp,
          col = "black", 
          size = 1,
           fill = "transparent") +
            coord_sf(xlim = c(metro_region_shp$centroid[[3]][1] - padding_width,
                              metro_region_shp$centroid[[3]][1] + padding_width - .5),
                     ylim = c(metro_region_shp$centroid[[3]][2] - padding_height,
                              metro_region_shp$centroid[[3]][2] + padding_height),
                    expand = FALSE)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

sgto_map2020 <- last_plot()

# 2022

col_pal <- c(#"#fef0d9",
             #"#fdcc8a",
             #"#fc8d59",
             "#e34a33",
             "#b30000"
             )

mob_indicators_1 %>% filter(REGION == "Metropolitana de Santiago" |
                                         REGION == "Valparaíso" | 
                                         REGION == "Libertador General Bernardo O'Higgins") %>% 
  filter(year == "2022") %>% 
  filter(distance_class == "<100") %>% 
   ggplot(aes(fill = new_jenk_class)) +
     geom_sf(col = "white", size = .1) +
   coord_sf() +
   scale_fill_manual(values= col_pal) +
  # scale_fill_brewer(palette = "OrRd", 
  #                   direction = 1, 
  #                   labels = c("10-421", "422-1231", "1,232-3,370", "3,371-22,761", "> 22,761")) +
   theme_map() +
   theme(plot.title = element_text(size = 60),
         legend.position = "none",
         panel.background=element_rect(colour="black")
         ) +
   labs(title = "Santiago") +
    geom_sf(data = metro_region_shp,
          col = "black", 
          size = 1,
           fill = "transparent") +
            coord_sf(xlim = c(metro_region_shp$centroid[[3]][1] - padding_width,
                              metro_region_shp$centroid[[3]][1] + padding_width - .5),
                     ylim = c(metro_region_shp$centroid[[3]][2] - padding_height,
                              metro_region_shp$centroid[[3]][2] + padding_height),
                    expand = FALSE)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

sgto_map2022 <- last_plot()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/chile/cloropleth-sgto2020.png", units="in", width=8, height=10, res=300)
  sgto_map2020
dev.off()

null device 
          1 

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/chile/cloropleth-sgto2022.png", units="in", width=8, height=10, res=300)

  sgto_map2022
dev.off()

null device 
          1 

Create bar plot

Short-distance

col_pal <- c(#"#fef0d9",
             #"#fdcc8a",
             "#fc8d59",
             "#e34a33",
             "#b30000"
             )

mob_indicators_1$end_polygon_name[mob_indicators_1$end_polygon_name == "San Felipe de Aconcagua"] <- "San Felipe"

# 2020
mob_indicators_1 %>% filter(year == "2020") %>%
  filter(distance_class == "<100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(end_polygon_name, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
       # axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_short2020 <- last_plot()

col_pal <- c(#"#fef0d9",
             #"#fdcc8a",
             #"#fc8d59",
             "#e34a33",
             "#b30000"
             )
# 2022
mob_indicators_1 %>% filter(year == "2022") %>%
  filter(distance_class == "<100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  distinct(end_polygon_name, .keep_all=TRUE) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(end_polygon_name, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
        #axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_short2022 <- last_plot()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/chile/barp_short2020.png", units="in", width=8, height=10, res=300)
  barp_short2020
dev.off()

null device 
          1 

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/chile/barp_short2022.png", units="in", width=8, height=10, res=300)

  barp_short2022
dev.off()

null device 
          1 

Long-distance

col_pal <- c("#fef0d9",
             "#fdcc8a",
             "#fc8d59"
             #"#e34a33",
             #"#b30000"
             )

# 2020
mob_indicators_1 %>% filter(year == "2020") %>%
  filter(distance_class == ">100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  distinct(end_polygon_name, .keep_all=TRUE) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(end_polygon_name, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
       # axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_long2020 <- last_plot()

# 2022

col_pal <- c(#"#fef0d9",
             #"#fdcc8a",
             "#fc8d59",
             "#e34a33"
             #"#b30000"
             )

mob_indicators_1 %>% filter(year == "2022") %>%
  filter(distance_class == ">100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  distinct(end_polygon_name, .keep_all=TRUE) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(end_polygon_name, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
       # axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_long2022 <- last_plot()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/chile/barp_long2020.png", units="in", width=8, height=10, res=300)
  barp_long2020
dev.off()

null device 
          1 

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/chile/barp_long2022.png", units="in", width=8, height=10, res=300)

  barp_long2022
dev.off()

null device 
          1 

Mexico

mobility data

# 2020 out
df20_out <- readRDS("/Volumes/RECAST/data/outputs/mexico/movements/2020_04_mov.rds") %>% 
  mutate(GEOMETRY = NULL) %>% 
  dplyr::filter(country == "MX") %>% 
  st_as_sf(coords = c("end_lon", "end_lat"), 
                                      crs = 'EPSG:4326')

# 2022 out
df22_out <- readRDS("/Volumes/RECAST/data/outputs/mexico/movements/2022_03_mov.rds") %>% 
  mutate(GEOMETRY = NULL) %>% 
  dplyr::filter(country == "MX") %>% 
  st_as_sf(coords = c("end_lon", "end_lat"), 
                                      crs = 'EPSG:4326')

boundary data

# admin boundaries shape
adm_shp <- st_read("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_MEX_shp/gadm41_MEX_2.shp")%>% 
  st_simplify(preserveTopology = T,
              dTolerance = 1000) %>% 
  st_make_valid() %>% 
  st_transform(crs = 'EPSG:4326')

Reading layer `gadm41_MEX_2' from data source 
  `/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_MEX_shp/gadm41_MEX_2.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 2457 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -118.3665 ymin: 14.53507 xmax: -86.71074 ymax: 32.71863
Geodetic CRS:  WGS 84

adm1_shp <- st_read("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_MEX_shp/gadm41_MEX_1.shp")%>% 
  st_simplify(preserveTopology = T,
              dTolerance = 1000) %>% 
  st_make_valid() %>% 
  st_transform(crs = 'EPSG:4326')

Reading layer `gadm41_MEX_1' from data source 
  `/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_MEX_shp/gadm41_MEX_1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 32 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -118.3665 ymin: 14.53507 xmax: -86.71074 ymax: 32.71863
Geodetic CRS:  WGS 84

# region shape
region_shp <- st_read("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_MEX_shp/gadm41_MEX_0.shp")%>% 
  st_simplify(preserveTopology = T,
              dTolerance = 1000) %>% 
  st_make_valid() %>% 
  st_transform(crs = 'EPSG:4326')

Reading layer `gadm41_MEX_0' from data source 
  `/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/data/maps/shapefiles/gadm41_MEX_shp/gadm41_MEX_0.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -118.3665 ymin: 14.53507 xmax: -86.71074 ymax: 32.71863
Geodetic CRS:  WGS 84

Data wrangling

# filter for start location as mexico city
df20_out <- df20_out %>%
  filter(start_polygon_name == 'Ciudad De México')

df22_out <- df22_out %>%
  filter(start_polygon_name == 'Ciudad De México')

Define distance

df20_out <- df20_out %>% mutate(
  distance_class = case_when(length_km < 100 ~ "<100",
                             length_km >= 100 ~ ">100"))

df22_out <- df22_out %>% mutate(
  distance_class = case_when(length_km < 100 ~ "<100",
                             length_km >= 100 ~ ">100"))

Sum of flows

# sum of outflows
outflows_df_20 <- df20_out %>% 
  filter(start_polygon_name != end_polygon_name) %>% 
  group_by(geometry, distance_class) %>% 
  dplyr::summarise(
    sum_outflow = sum(n_crisis, na.rm = T)) %>% 
  ungroup() 

`summarise()` has grouped output by 'geometry'. You can override using the `.groups` argument.

outflows_df_22 <- df22_out %>% 
  filter(start_polygon_name != end_polygon_name) %>% 
  group_by(geometry, distance_class) %>% 
  dplyr::summarise(
    sum_outflow = sum(n_crisis, na.rm = T)) %>% 
  ungroup() 

`summarise()` has grouped output by 'geometry'. You can override using the `.groups` argument.

Check alignment between data frames

# plot boundaries - test 
p <- ggplot() + 
  geom_sf(data = adm_shp,
          color = "gray60", 
          size = 0.1) 

last_plot()

# check geometry fits into either admin boundaries or region boundaries
p <- p +
  geom_point(data = outflows_df_20,
    aes(geometry = geometry),
    stat = "sf_coordinates"
  ) 
last_plot()

Add year

# assign a year to each
outflows_df_20$year <- '2020'
outflows_df_22$year <- '2022'

Bind outflows

outflows_df <- rbind(outflows_df_20, outflows_df_22)

Join flows to polygons

mob_indicators_1 <- st_join(adm_shp, outflows_df)

classify data into quantiles

mob_indicators_1$new_jenk_class <- classify_intervals(mob_indicators_1$sum_outflow, n = 5, style = "quantile", factor = TRUE)

Warning: var has missing values, omitted in finding classes

outflow_labels_1 <- levels(mob_indicators_1$new_jenk_class)
outflow_labels_1 <- gsub("^.|.$", "",  outflow_labels_1)
outflow_labels_1 <- gsub("\\.[0-9]+", "",  outflow_labels_1)
outflow_labels_1 <- gsub("\\,", "-",  outflow_labels_1)
levels(mob_indicators_1$new_jenk_class) <- outflow_labels_1

# mob_indicators_1 <- mob_indicators_1 %>% 
  # mutate(new_jenk_class = str_replace(new_jenk_class, ",", "-"))

# change geometry
shp_reg <- region_shp %>% st_transform(crs = 'EPSG:4326')

mob_indicators_1$new_jenk_class <- mob_indicators_1$new_jenk_class %>% replace_na("10-92")

mob_indicators_1 <- na.omit(mob_indicators_1)

Plotting

Creating individual country map

shp_reg$centroid <- shp_reg %>% 
  st_centroid() %>% 
  st_geometry()

Warning: st_centroid assumes attributes are constant over geometries

padding_width <- 17 
padding_height <- 22

new_y_labels <- c("Short-distance (<100km)", "Long-distance (>100km)")
names(new_y_labels) <- c("<100", ">100")

ggplot(data = mob_indicators_1, 
                           aes(fill = new_jenk_class )) +
  geom_sf(col = "transparent", size = .1) + 
  coord_sf() +
  scale_fill_brewer(palette = "OrRd", 
                    direction = 1, 
                    labels = c("10-92", "93-443", "444-1,654", "1,655-4,458", "> 4,458"), 
                    na.value="black") +
  facet_grid(distance_class ~ year, 
             labeller = labeller(distance_class = new_y_labels)) +
  labs(title = "C. Mexico",
       fill = "Number of out-moves") +
  theme_map() +
  guides(
    color='none',
    fill = guide_legend(
      keywidth = 4, 
      keyheight = 1,
      nrow = 1,
      title.position="top",
      label.position="bottom"
    )
  ) +
  theme(plot.title = element_text(size = 16, face = "bold"),
        plot.margin = margin(1,0,1,0,"cm"),
        legend.title = element_markdown(
          size=10, 
          face = "bold", 
          hjust=0.5, 
          lineheight=0.45,
          color="black",
          margin=margin(0,0,-0.2,0,"cm")
          ),
        legend.text = element_text(size = 8),
        legend.position = "bottom",
        legend.spacing.x = unit(0, 'cm'),
        panel.background = element_rect(fill = "gray99", colour = "gray99")
        ) +
    geom_sf(data = shp_reg,
          col = "black", 
          size = 1,
          fill = "transparent") +
            coord_sf(xlim = c(shp_reg$centroid[[1]][1] - padding_width, 
                              shp_reg$centroid[[1]][1] + padding_width), 
                     ylim = c(shp_reg$centroid[[1]][2] - padding_height - 7, 
                              shp_reg$centroid[[1]][2] + padding_height - 2), 
                    expand = FALSE)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

outflow_plot_mex <- last_plot()

Saving country map

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/mexico/cloropleth-mex1.png", units="in", width=8, height=10, res=300)
  outflow_plot_mex
dev.off()

null device 
          1 

Create map for Mexico City

Short-distance

mxcity_region_shp <- adm1_shp %>% filter(NAME_1 == "Distrito Federal" |
                                           NAME_1 == "México" |
                                           NAME_1 == "Hidalgo" |
                                           NAME_1 == "Puebla" |
                                           NAME_1 == "Morelos" |
                                           NAME_1 == "Guerrero")

mxcity_region_shp$centroid <- mxcity_region_shp  %>% 
  st_centroid() %>% 
  st_geometry()

Warning: st_centroid assumes attributes are constant over geometries

padding_width <- 4.5
padding_height <- 3.5

col_pal <- c("#fef0d9",
             #"#fdcc8a",
             #"#fc8d59",
             #"#e34a33",
             "#b30000")

# 2020
mob_indicators_1 %>% filter(NAME_1 == "Guerrero" |
                              NAME_1 == "Hidalgo" |
                              NAME_1 == "México" |
                              NAME_1 == "Puebla" |
                              NAME_1 == "Tlaxcala") %>% 
  filter(year == "2020") %>% 
  filter(distance_class == "<100") %>% 
   ggplot(aes(fill = new_jenk_class)) +
     geom_sf(col = "white", size = .1) +
   coord_sf() +
   scale_fill_manual(values= col_pal) +
  # scale_fill_brewer(palette = "OrRd", 
  #                   direction = 1, 
  #                   labels = c("10-421", "422-1231", "1,232-3,370", "3,371-22,761", "> 22,761")) +
   theme_map() +
   theme(plot.title = element_text(size = 60),
         legend.position = "none",
         panel.background=element_rect(colour="black")
         ) +
   labs(title = "Mexico City") +
    geom_sf(data = mxcity_region_shp,
          col = "black", 
          size = 1,
           fill = "transparent") +
            coord_sf(xlim = c(mxcity_region_shp$centroid[[1]][1] - padding_width,
                              mxcity_region_shp$centroid[[1]][1] + padding_width - .5),
                     ylim = c(mxcity_region_shp$centroid[[1]][2] - padding_height,
                              mxcity_region_shp$centroid[[1]][2] + padding_height),
                    expand = FALSE)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

mxcity_map2020 <- last_plot()

# 2022

col_pal <- c("#fef0d9",
             "#fdcc8a",
             "#fc8d59",
             "#e34a33",
             "#b30000")

mob_indicators_1 %>% filter(NAME_1 == "Guerrero" |
                              NAME_1 == "Hidalgo" |
                              NAME_1 == "México" |
                              NAME_1 == "Puebla" |
                              NAME_1 == "Morelos" |
                              NAME_1 == "Tlaxcala") %>% 
  filter(year == "2022") %>% 
  filter(distance_class == "<100") %>% 
   ggplot(aes(fill = new_jenk_class)) +
     geom_sf(col = "white", size = .1) +
   coord_sf() +
   scale_fill_manual(values= col_pal) +
  # scale_fill_brewer(palette = "OrRd", 
  #                   direction = 1, 
  #                   labels = c("10-421", "422-1231", "1,232-3,370", "3,371-22,761", "> 22,761")) +
   theme_map() +
   theme(plot.title = element_text(size = 60),
         legend.position = "none",
         panel.background=element_rect(colour="black")
         ) +
   labs(title = "Mexico City") +
    geom_sf(data = mxcity_region_shp,
          col = "black", 
          size = 1,
           fill = "transparent") +
            coord_sf(xlim = c(mxcity_region_shp$centroid[[1]][1] - padding_width,
                              mxcity_region_shp$centroid[[1]][1] + padding_width - .5),
                     ylim = c(mxcity_region_shp$centroid[[1]][2] - padding_height,
                              mxcity_region_shp$centroid[[1]][2] + padding_height),
                    expand = FALSE)

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

mxcity_map2022 <- last_plot()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/mexico/cloropleth-mxcity2020.png", units="in", width=8, height=10, res=300)
  mxcity_map2020
dev.off()

null device 
          1 

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/mexico/cloropleth-mxcity2022.png", units="in", width=8, height=10, res=300)

  mxcity_map2022
dev.off()

null device 
          1 

Create bar plot

Short-distance

# 2020
col_pal <- c("#fef0d9",
             #"#fdcc8a",
             #"#fc8d59",
             #"#e34a33",
             "#b30000")

mob_indicators_1 %>% filter(year == "2020") %>%
  filter(distance_class == "<100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(NAME_2, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
       # axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_short2020 <- last_plot()

# 2022

 col_pal <- c(#"#fef0d9",
#              "#fdcc8a",
#              "#fc8d59",
#              "#e34a33",
             "#b30000")

mob_indicators_1 %>% filter(year == "2022") %>%
  filter(distance_class == "<100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  distinct(NAME_2, .keep_all=TRUE) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(NAME_2, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
        #axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_short2022 <- last_plot()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/mexico/barp_short2020.png", units="in", width=8, height=10, res=300)
  barp_short2020
dev.off()

null device 
          1 

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/mexico/barp_short2022.png", units="in", width=8, height=10, res=300)

  barp_short2022
dev.off()

null device 
          1 

Long-distance

# 2020
col_pal <- c(#"#fef0d9",
             #"#fdcc8a",
             #"#fc8d59",
             "#e34a33",
             "#b30000")

mob_indicators_1 %>% filter(year == "2020") %>%
  filter(distance_class == ">100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(NAME_2, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
       # axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_long2020 <- last_plot()

# 2022
col_pal <- c(#"#fef0d9",
             #"#fdcc8a",
             #"#fc8d59",
             #"#e34a33",
             "#b30000")

mob_indicators_1 %>% filter(year == "2022") %>%
  filter(distance_class == ">100") %>% 
  arrange(desc(sum_outflow)) %>% 
    mutate(outflow_total = sum(sum_outflow),
         outflow_percent = round( (sum_outflow/outflow_total)*100, 2 ) ) %>% 
  head(n = 10) %>% 
  ggplot(aes(x = reorder(NAME_2, outflow_percent), y = outflow_percent, fill = new_jenk_class)) + 
  geom_bar(stat="identity" ) + 
  scale_fill_manual(values= col_pal) +
  theme_tufte2() +
  coord_flip() +
  labs(title = "Top ten outflows",
       y = "Percent (%)") +
  theme(plot.margin=margin(1,0,1,1,"cm"),
        legend.position = "none",
        plot.title = element_text(size = 35),
        axis.text = element_text(size = 32),
        axis.title.y = element_blank(),
       # axis.title.x = element_blank(),
        axis.ticks.y = element_blank()) # Set custom y-axis labels

barp_long2022 <- last_plot()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/mexico/barp_long2020.png", units="in", width=8, height=10, res=300)
  barp_long2020
dev.off()

null device 
          1 

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/mexico/barp_long2022.png", units="in", width=8, height=10, res=300)

  barp_long2022
dev.off()

null device 
          1 

Saving combined maps

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/cloropleth-map.png", units="in", width=8, height=10, res=300)
  outflow_plot_arg + outflow_plot_chi + outflow_plot_mex
dev.off()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/cloropleth-map1.png", units="in", width=8, height=10, res=300)
  wrap_plots(outflow_plot_arg,
             outflow_plot_chi,
             outflow_plot_mex)
dev.off()

png("/Users/franciscorowe/Dropbox/Francisco/Research/in_progress/recast/cepal-report/outputs/cloropleth-maps/cloropleth-map2.png", units="in", width=8, height=10, res=300)
  plot_grid(outflow_plot_arg, outflow_plot_chi, outflow_plot_mex, ncol=3, align="h")
dev.off()